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In the replica-exchange molecular dynamics method, where constant-temperature molecular dy- 
namics simulations are performed in each replica, one usually rescales the momentum of each par- 
£N| ' ticle after replica exchange. This rescaling method had previously been worked out only for the 

» | | Gaussian constraint method. In this letter, we present momentum rescaling formulae for four other 

f"S |. commonly used constant-temperature algorithms, namely, Langevin dynamics, Andersen algorithm, 

1 Nose-Hoover thermostat, and Nose-Poincare thermostat. The effectiveness of these rescaling meth- 

ods is tested with a small biomolecular system, and it is shown that proper momentum rescaling is 
s ! , necessary to obtain correct results in the canonical ensemble. 

PACS numbers: 

^ . Monte Carlo (MC) and molecular dynamics (MD) simulations with generalized-ensemble algorithms have been 
widely used for studies of proteins and peptides (for a review, see, e.g., ref. 1). Among generalized-ensemble algorithms, 
the replica-exchange method (REM) Q (the method is also referred to as parallel tempering Q) is a useful simulation 
method because there is no need to determine the weight factor before the simulation. The original REM was proposed 
for MC simulations and momenta of particles do not have to be considered, but momenta have to be included in 
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replica-exchange molecular dynamics (REMD) simulations and are usually rescaled when the exchange is performed. Q 
Although only dynamical variables of a physical system are considered and the usual Boltzmann weight factor is used 
in the detailed balance condition in the original REMD,[i[ the states of the system can be specified by coordinates 
and momenta of atoms and additional dynamical variables, or the distribution function can be different from that of 
the canonical ensemble in some constant temperature algorithms. For example, the Nose-Hoover thermostat has a 
velocity of the thermostat, and the distribution function includes a kinetic energy term of this velocity. [f| Thus, one 
has to use a different rescaling method in REMD simulations depending on different constant temperature algorithms. 

Let us consider a REM system, which consists of M non-interacting replicas of the original system in the canonical 
ensemble at M different temperature values T m (m = 1, • ■ • , M), and let i (i = 1, ■ • • , M) be a label which specifies 
7— I ■ the replica. Then there is a one-to-one correspondence between the replica label and the temperature label and we 
^ ' can introduce a permutation function a and the inverse function er" 1 as 

t> : r 

IT), \i = a(m), m = <7 

^ '. \j = (r(n), n = a- 1 (j), 

. where i and j stand for the replica labels and m and n stand for the temperature labels, respectively. The states 
of replica i at temperature value T m are specified by Xm = (q^,p^,a^) m . Here, = {q± ,q^j} and = 

{Pi\'" )Pjv) are coordinates and momenta of N atoms in replica i, respectively, and is a set of additional 
• • , dynamical variables depending on constant temperature algorithms in replica i. The states of the REM system are 
.> specified by X = {x [ ? W \ - ■ ■ ,x$ M)] }. 

Suppose a pair of replicas i and j are exchanged and let X' be a state after the replicas are exchanged. Namely, 
5h ' we try to change the state of the REM system as follows: 

X = {■ ■ ■ , x m , • • ■ , x^ ,•••}—>■ X = {■••, x^ , ■ ■ ■ , x n , ■ • ■ }, (2) 



where Xm' and in' are defined by 



'x#' = (g&],p&]', a y]') 



(3) 



respectively, and explicit forms of the variables and a^']' are determined below. After this transitions 

of states, the permutation function a is replaced by a new permutation function a' defined by 

\i = a'(n). 1 ^ 
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To make the REM system approach an equilibrium state, the detailed balance condition is imposed: 

P (X) w (X -> X') = P (X') w (X' -> X) , 
where w (X — > X') is the transition probability from X to X' and P (X) is given by 

M 



(5) 



(6) 



m— 1 



Here, f m is a distribution function at temperature T m . The transition probability w(X —> X') is then obtained by 
the Metropolis criterion: [6| 



w(X -4 X') 



1, 



P(X) 



(7) 



In general, the transition probability will be different from that of the replica-exchange MC algorithm, which is 
given by Q 



■ (X -> X') = min [1, exp(-zA)] , 



where 



4 = (A - Pm) 



E(q [ 



E{qW) 



(8) 



(9) 



/3 m is the inverse temperature defined by f3 m = \jk^T m (ks is the Boltzmann constant), and E is the potential 
energy function. However, we want this probabilty w(X — > X') in eq. of replica exchange to be independent 
of the constant-temperature algorithms. One of such a choice is to use also the same probability as in eqs. © 
and ([9]). We believe that this is the most natural choice. Therefore, in the following discussion, we determine how 
to rescale the variables and a^'' for various constant temperature algorithms so that the transition 

probability is given by eqs. ((U and ((9]). The constant temperature algorithms that are discussed in this letter are 
the Gaussian constraint method. B-jl Langevin dynamics, Andersen algorithm, [llj Nose-Hoover thermostat, Q 
and Nose-Poincare thermostat. |12| 

The Gaussian constraint method, in which the states are specified by x = (q,p), has the following distribution 
function: Q 



f{q,p) oc S (j2 |L - exp [-0E(q)] , 



(10) 



where g — 3N — 1 if there are no constraints in the system. Although the distribution function is not that of the 
canonical ensemble because the kinetic energy is fixed, it can be shown that the proper rescaling method in eqs. ((5]) 
and © is given by [3] 



,[*]' 




21ml*] 




(11) 



This is the original momentum scaling introduced in the REMD method. Q| 

In the Langevin dynamics and the Andersen algorithm, the states are also specified by x = (q,p) and the distribution 
function is that of the canonical ensemble, [13, [ll| that is 



f(q,p) oc exp 



JY 



Pi 

2mi 



E(q) 



(12) 



The rescaling method in the algorithms is obtained, following the original REMD paper Q and again given by eq. 



CD 



In the Nose-Hoover thermostat, the states are specified by x = (q,p, C) and the distribution function is given by 5] 



f{q,P,0 oc exp 



(13) 
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where £ is a thermostat velocity and Q is its mass parameter. The mass parameter can have different values in each 
replica in REMD simulations. It is found that the transition probability is given by eqs. © and ©, if the momenta 
are rescaled by eq. (fTT|) and C^' and are rescaled by 



y T^vaQn y ^nQm 

where Q m and Q n are the mass parameters in the replicas at temperature values T m and T„, respectively. The 
rescaling method can be generalized to the Nose-Hoover chains [H| in a similar way. 

In the Nose-Poincare thermostat, the states are specified by x = (q,p,s,ir) and the distribution function is given 

by 

f[q,P, s,tt) oc S [s(H n -£)] , (15) 
where £ is an initial value of and is the Nose Hamiltonian, [3, [H[ which is given by 

N ~2 2 

^ = E^2+^) + ^+. 9 fcBTlog, (16) 

Here, g(= 3N) is the number of degrees of freedom, s is a position variable of the thermostat, it is a momentum 
conjugate to s, and pi is a virtual momentum, which is related to the real momenta pt as pi — pi/s. A rescaling 
method of the Nose-Poincare thermostat can be given by eq. (fTTj) and 



^' = J^r^ ] , ^ ] ' = \ ^r^\ (17) 



s w = s w exp 



g\j]> — s [j] g X p 



J_ / E(q®)-£„ 

gkB \ 



E(qW)-£ n 

T T 

J_ / E(qW)-£ n E{qW)-£ n 
gk B V T n 



T 



(18) 



where £ m and £ n are initial values of An in the simulations with T m and T n , respectively. Note that the real momenta 
have to be used in the rescaling method in eq. (|11[) . not the virtual momenta. 

In order to test the validity of the rescaling methods, we performed REMD simulations with the five constant tem- 
perature algorithms. We used a system of a short peptide, Met-enkephalin, in gas phase, which has the amino-acid 
sequence Tyr-Gly-Gly-Phe-Met. The N-terminus and the C-terminus of the peptide were blocked by the acetyl group 
and the N-methylamide group, respectively, and the initial structure was a fully extended one in all the replicas. We 
used the following eight temperature values: (T 1( --- ,T 8 ) = (200,239,286,342,409,489,585,700) in kelvin. Initial 
velocity of each atom was given by the Maxwell-Boltzmann distribution corresponding to the temperature in each 
replica. The simulations were performed by the tinker program package. [16| Several of the programs in the pack- 
age were modified and a few programs were added so that REMD simulations with the five constant temperature 
algorithms can be performed. We used the AMBER parm99 force field pj]], and the dielectric constant was set to 
1.0. The unit time step was set to 0.5 fs and the simulation time was set to 5.0 ns for each replica. We tried to 
exchange the replicas at every 10 fs and the data of the simulations were sampled just before the exchange trials. 
The replica exchange was performed by exchanging the temperatures instead of the dynamical variables. The pairs 
of replicas corresponding to neighboring temperature were simultaneously exchanged and the two choices of pair- 
ing, (Ti, T2), (T3, T4), • • • and (T2, T3), (T4, T5), • • • , were used alternately. These conditions were common to all the 
simulations. The computational details for each constant temperature algorithm are described below. 

In the Gaussian constraint method, the integration method proposed by Zhang [18j was used. In the Langevin 
dynamics we used the integration method proposed by Mannella [19| and the friction coefficient was set to 5.0 ps _1 . 
In the Andersen algorithm the Velocity Verlet algorithm was used and the random collision frequency was set to 0.1 
fs _1 -atom _1 . In the Nose-Hoover thermostat we used the integration method proposed by Martyna et al. [20[ and 
the values of mass parameters Q were determined by the following equation [15j : 



(19) 
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FIG. 1: Potential energy distributions in the simulations with and without REM for the Nose-Poincare thermostat. The results 
with REM are represented by the solid curves and without REM by the crosses. The distributions of the right side correspond 
to the highest temperature (700 K) and those of the left side correspond to the lowest temperature (200 K). 

where oj = 2it/t is a frequency of the thermostat and we set r to 0.01 ps. The initial values of £ was set to 0.0. In the 
Nose-Poincare thermostat we used the symplectic integrator proposed by Nose [2l[ , which was recently generalized to 
biomolecular systems which include rigid-body molecules. [22[ The value of each mass parameter was determined by 
the same way and had the same value as the Nose-Hoover thermostat. The initial values of s and ir were set to 1.0 
and 0.0, respectively. For practical purposes, eq. (TT5)) cannot be used directly because the Nose Hamiltonian is not 
actually conserved in the simulations. Thus, modifications of eq. (fT5)) are required. A simple solution is to replace 
£ m for replica i and £ n for replica j by the current value of the Nose Hamiltonian. This prescription works well as is 
shown below. 

We also performed conventional canonical simulations (that is, without REM) with five constant temperature 
algorithms under the same conditions as in the REMD simulations so that the results with REM may be compared 
to those without REM. Figure [T] shows the distributions of the potential energy in the simulations with and without 
REM for the Nose-Poincare thermostat. At the highest temperature the distribution with REM agrees with that 
without REM, while at the lowest temperature the potential energy of the simulation with REM was lower than 
without REM. Essentially the same results were also obtained for the other constant temperature algorithms and 
the obtained distributions for all the REMD simulations agreed with each other, which are shown in Fig. [2j The 
results at the lowest temperature by the canonical simulation without REM in Fig. [T]thus are wrong because it got 
trapped in states of local-minimum energy. These results imply that REMD simulations with the rescaling methods 
can generate the correct canonical distributions. Table U lists the acceptance ratios for each REMD simulation. The 
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FIG. 2: (Color online) Potential energy distributions that were obtained by the REMD simulations with the five constant 
temperature algorithms. All of the results are plotted and essentially coincident among the five algorithms. 
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TABLE I: Acceptance ratios of replica exchange between pairs of the temperature values in the REMD simulations for all the 
constant temperature algorithms. In this table, G, L, A, NH, and NP stand for Gauss, Langevin, Andersen, Nose-Hoover, and 
Nose-Poincare, respectively. 





(Ti,T a ) 


(T 2 ,T 3 ) (T 3 ,T 4 ) 


(24, T 5 ) 


(T 5 ,T 6 ) (T 6 ,T 7 ) (T 7 ,T 8 ) 


G 


0.148 


0.140 


0.137 


0.132 


0.129 


0.131 


0.139 


L 


0.149 


0.142 


0.142 


0.128 


0.128 


0.129 


0.134 


A 


0.151 


0.143 


0.141 


0.130 


0.125 


0.128 


0.135 


NH 


0.151 


0.144 


0.139 


0.130 


0.127 


0.132 


0.138 


NP 


0.151 


0.144 


0.142 


0.132 


0.127 


0.128 


0.136 



acceptance ratios are almost the same for all the simulations and therefore it is found that the efficiency of REMD 
simulations is independent of the kind of constant temperature algorithms. Figure [3] shows the structures of the 
lowest potential energy state in all the REMD simulations. Although the orientations of some side-chains are slightly 
different among the structures, the backbone structures were obtained with almost the same accuracy (less than 0.6 
A in rms deviation). 

To show that the rescaling methods are necessary to perform proper REMD simulations, we performed a REMD 
simulation with the Nose-Poincare thermostat, in which eqs. ((5J) and ^ were used but the dynamical variable s was 
not rescaled. Figure 0] shows the time series of the potential energy in the simulations when the dynamical variable 
s was rescaled and not rescaled. If s was rescaled properly, the REMD simulation worked properly but otherwise the 
potential energy tended to diverge. From these results, it is found that if dynamical variables are not properly rescaled 
in REMD simulations, then the simulation sometimes cannot provide correct results. Therefore, one has to use the 
appropriate rescaling methods for each constant temperature algorithm when REMD simulations are performed. 

In this letter, we proposed the rescaling methods which treat dynamical variables properly in REMD simulations for 
various constant temperature algorithms. REMD simulations with the proper rescaling methods can provide correct 
canonical-ensemble distributions. With these rescaling methods for familiar constant temperature algorithms, REMD 
methods will become more applicable for simulations of various molecular systems. 
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FIG. 3: (Color online) Structures of the lowest potential energy state that were obtained in the REMD simulations with (a) 
the Gaussian constraint method, (b) the Langevin dynamics, (c) the Andersen algorithm, (d) the Nose-Hoover thermostat, and 
(e) the Nose-Poincare thermostat, respectively. These figures were created by the vmd software. [23| 
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FIG. 4: (Color online) Time series of the potential energy in the REMD simulations at 700 K with the Nose-Poincare thermostat. 
The results with the rescaling method are stable and flat (red curve), and those without the rescaling method are unstable 
(green curve). 



